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Abstract 

The coupled strongly implicit method described in Part 1 of 
this study is combined with a deferred-corrector spline solver 
for the vortici ty-stream function form of the Navier-Stokes oqua 
tion. Solutions for cavity, channel and cylinder flows are ob- 
tained with the fourth-order spline 4 procedure. The strongly 
coupled spline corrector method converges as rapidly as the 
finite difference calculations of Part 1 and also allows for 
arbitrary large time increments for the Reynolds numbers con- 
sidered, (R^£l000) . In some cases fourth-order smoothing or 
filtering is required in order to suppress high frequency oscil- 
lations . 


1. INTRODUCTION 

In recent papers the present authors have formulated 
various higher-order collocation techniques based on polynomial 
interpolation or Hermitian discretization procedures. Three-point * 
formulations leading to fourth-order and sixt.i-order methods were 
considered. In order to obtain equal accuracy, these higher-or- 
der methods require fewer mesh points when compared with second- 
order accurate I i;iite-dif ference techniques. This can moan a re- 
duction in computer storage and time. In particular, the fourth- 
order accurate spline 4 method generally requires one-quarter 

*This "research was sponsored by the National Aeronautics and Space 
Administration, Langley Research Center Hampton, Va., under 
Grant NSG-1244 


1 


the number of points, in each coordinate direction, compared 
to finite-difference solutions of comparable accuracy. 

In addition to improvements in accuracy, these methods have 
certain desitable properties. The discrete equations remain 
block- tridiagonal in character. In vievv of the Hermitian formu- 
lation, where the derivatives are treated as unknowns, the appli- 
cation of derivative boundary conditions is somewhat more 
straightforward, and with non-uniform grids the decerioration in 
accuracy is less severe than with conventional finite-difference 
methods. However, the resulting system of equations is more 
complex than those of the usual second-order accurate finite-dif- 
ference procedures. For a single differential equation, finite- 
differences lead to a scalar tridiagbnal system, whereas the 
higher order techniques are block tridiagonal; specifically, 2x2 
blocks are required for the fourth order spline metliods and 3x3 
blocks for the sixth order Hermite systems. Consequently, for a 
given mesh the storage requirements become larger. The solution 
procedure therefore requires more operations and is more time 
consuming than the scalar algorithm usually applied for the tri- 
diagonal finite-difference equations. However, for equal accur- 
acy the spline methods require fewer points and therefore are 
still more efficient. 

( 4 ) 

In a recent study ' some of the complexities associated 
with computer storage and the general spline solution procedure 
have been eliminated. This was accomplished by "uncoupling" the 
higher-order correction from the coupled lower-order ("finite- 
difference") system of equations. This deferred-corrector method 
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can therefore result in a significant reduction in computer 

storage. A 2x2 block tridiagonal system is reduced to two scalar 

tridiagonal systems; 3x3 blocks reduce to a tridiagonal and a 

2x2 block tridiagonal system. In particular, the method can be 

programmed into a subroutine package to be used with existing 

finite-difference codes in order to render higher-order accuracy 

to the final solutions. This simplicity is achieved by treating 

( 

the higher-order accurate terms as explicit correctors to a 

modified finite-difference equation. This correction is evaluated 

in such a manner that the resulting procedure is consistent and 

( 4 ) 

unconditionally stable. Related non-iterative methods of de- 

ferred-correction have been discussed by several authors fer or- 
dinary differential equations. If these procedures are ex- 
tended so as to achieve full convergence, for ODE's or for PDE's 
in a temporal or spatial marching procedure, an instability can 
result. The analysis of reference (4) corrects this deficiency. 
Finally, it is noteworthy that the Khosla-Rubin difference meth- 
od, which provides an explicit central-difference corrector to 
an implicit upwind-difference scheme, is in fact a lower-order 
unconditionally stable adaptation of the deferred-corrector method. 

It should be noted that by uncoupling the spline or Hermite 
corrector and thereby reducing the size of the matrix blocks, it 
is now possible to obtain higher-order solutions and still main- 
tain a considerable degree of coupling between the dependent 
variables in the governing systems of equations. In effect, the 
equations arc coupled solely through the implicit or finite-dif- 
fercnce-like terms. This is particularly important for the Navier- 
Stokes calculations, as the rate of convergence is increased 
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notably when the equations are coupled. The rate of convergence 

docs not ippoar to be as sensitive to the uncoupling of the 

higher-order spline terms. In this way wo have been able to 

develop a coupled solver for the vorticity (oj) -stream function 

(ij;) description of the Navier-Stokes equations. Only (2x2) 

blocks must be inverted in order to obtain fourth-order accuracy. 

This procedure has been used with ADI, predictor-corrector, direct 

solver and strongly implicit two-dimensional solution tech- 
(4 ) 

' niques. In the present study the coupled strongly implicit 

method is used exclusively to treat a variety of problems. 

In Part 1^^^ of the present study the coupled strongly im- 
plicit method has been detailed. The importance of implicitly 
coupling the boundary conditions has been demonstrated. Second- 
order accurate central-finite-dif ference solutions have been ob- 
tained for (i) the flow in a driven cavity, (ii) the temperature 
and heat transfer in the cavity, (iii) a channel wi'th a rearward 
facing step, and (iv) a circular cylinder with a splitter plate. 

In the second part of the analysis presented here, these geome- 
tries are reconsidered with the deferred-corrector spline adaption 
of the coupled strongly implicit method. 

In section 2 the spline collocation equations and analysis are 
reviewed. In section 3 the deferred-corrector technique of ref- 
erence 4 is described for both the K-R and spline formulations. 

In section 4 the coupled strongly implicit spline corrector mJithod 
is reviewed for the i|/-o) Navier-Stokes equations. Finally, the 
spline solutions for the problems of Part 1 are presented. 
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2. Spline 4 Formulation 

A variety of fourth and sixth-order spline and hermite pro- 
cedures have been derived in references 1 and 2. The major ad- 
vantages of these techniques over conventional three-point 
finite-difference methods arc also discussed in these papers. 

In the present study we will be concerned only with the spline 4 
procedure. The following two- and three-point spline relation- 
ships between the variable u^ and its first (t^) and second 

derivatives (L^) in the y-coordinate direction (i fixed) are re- 

(1 2 ) 

quired for the execution of the iterative procedure. ' 
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where o =* h.,,/h. and h. is the mesh width. 
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The truncation errors of the derivatives are given as: 
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Similar expressions arc obtained for m. M. representing (u ) , (u ) , 

X 0 X X XX 

respectively . 

It should be noted that with o ■ 1+0 (h^), even the variable 
grid representation leads to fourth-order accuracy. For given 
values of u^, a single inversion of a scalar tridiagonal matrix 
(equation (Id)} is required for the evaluation of f ^ and . If 
Uj is to be evaluated from the solution of a differential equation, 
the spline 4 procedure leads to a 3x3 bloc)<-tridiagonal system 
for (u,f ,K)^; this can always be reduced to a 2x2 system with the 
help of equations (lb). In references 1 and 2 , this coupled pro- 
cedure was used in the solution of a variety of flow problems. In 
reference 4, the spline 4 method was further simplified such that 
the 2x2 system for (u,K)j is reduced to the inversion of two scalar 
tridiagonal systems for Uj,Kj. This reduces the computer storage 
considerably and allows the higher-order improvement to bo pro- 
grammed as a simple corrector to a finite-difference code. 

It is this deferred-corrector approach that is applied here 
for the solution of the (i|(-aj) Navier-Stokes equations. The pro- 
cedure is completely general and can be extended to any higher- 
order method. In reference 3 solutions are obtained with the 
Hormite 6 procedure for turbulent boundary layers. The stability 
analysis is reported in reference 4. 


6 


3. Iterative Solution Procedure: Equations 

A general deferred-corrector procedure for several higher- 

order numerical methods has been described in reference 4, where 

the necessary stability conditions have also been presented. In 

this study, we focus specifically on the spline 4 procedure, 

applied as a deferred-corrector, for the form of the Navier- 

Stokes equations. This approach is unconditionally stable, differs 

(5) 

from earlier methods of deferred-correction and represents a 
natural extension of the second order K-R scheme to higher-or- 
der differencing procedures. We are primarily concerned here with 
the steady-state solution procedure; however, the method is des- 
cribed for the general unsteady equations. For the applications 
described herein, large time steps (AslO®) are used in all of the 
calculations. 

For the Navier-Stokes equations, the general iterative 

procedure can be written, in the transformed variables (x,y) , as 
follows : 

Vorticity Transport Equation ; 

The difference approximation for the vorticity equation (3) 

j <D^ + ^“‘*’^x * ^ 

* e 
is 
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H ,m,£. ,ih» Lf M are the spline approximations to w ,u) , fvto) (uu) , u) 

y X y X yy 

and u) , respectively, see section 2. J is the Jacobian of the 

mapping function to the (x,y) plane. The transformations will 

bo reviewed for the various problems as they are presented, 

(7) 

see Part 1. The constant c is prescribed as zero or unity; 

A and B are constants which are chosen in order that the itera- 


(4a ) 


(4b ) 
(4c ) 

(4d) 

(4e) 

(4f,g) 
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tivc procedure is unconditionally stable. When complete con- 
vergence is achieved, the terms involving A and B cancel, so 
that with c ■ 1, a higher-order numerical solution of the 
Navier-Stokes equations is recovereo. 

The convective terms on the left-hand side of equation 
(4a) are of the form obtained with an upwind finite-difference 
discretization. For A»B“1 the left-hand side corresponds 
exactly with upwind differencing for convection and central dif- 
ferences for diffusion. The higher-order terms appear as ex- 
plicit correctors on the righthand side. If the spline correc- 
tors were replaced with orly second-oider correctors the K-R 
method is recovered. 

This implicit-explicit splitting of the convective terms 

is not unique, and is prescribed here to provide the appropriate 

five-point formulas required for the coupled 2x2 algorithm pre- 

(7) 

sented in Part I of this study. The novel feature of the 

present procedure is the introduction of r.he A, B splitting. 

This is not found in the usual deferred-corrector methods. 

The final A,B values are specified from stability considerations. 

(4 ) 

The stability cinalysis provides the minimum relative weights 
of implicit convection (A) and diffusion (B) required to achieve 
unconditional stability. For second-order schemes, equal weight- 
ing with A»=B=1 is unconditionally stable for i.e., multiple 

iteration is convergent with the K-R method. On the contrary, in 
order to achieve higher-order accuracy, A and B must be non-equal 

( 4 \ 

and different from unity. For the spline 4 method, A«*2 and 

B’»3 are the minimum integer values that load to an unconditionally 
t table corrector. 
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In certain eases it is desirable to use the deferred-cor- 
rector in a step-like fonliLon. An initial first-order accurate 
solution is obtained with A«D=1, r«0, i.e., upwind-differencing. 
This can be corrected to second-order accuracy with the K-R 
corrector. Finally, with C"l, by selecting appropriate values 
of A and D, the order of the solution can be further upgraded. 

For spline 4, A*2, D=3 and the spline formulas of section 2 are 

I 

applied. For some other higher-order methods the A,B values 
are given in reference 4. In this manner, a first-order accurate 
solution with considerable artificial viscosity is corrected to 
sccond-and then fourth-order. Alternatively, one could proceed 
directly from first to fourth-order or directly from second to 
higher-order, see reference 4. The higher-order solutions have 
little numerical diffusion. 

Stream Function Equation 

The stream function equation (5) 

V^i|/ - u)/J (i) 

is approxii.iated by 
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^ij and M... are the spline apt- roximat ions to and , respec- 
ij yy ^ XX ^ 

tivcly. The velocities u and v are related tr* and by the 

the spline relationships given in section 2; the constant C plays 

a role similar to that of B in (4a) and c is zero or one as be* 

fore. For finite-differences C ■ 1 and for spline 4,C^1. For 

( 4 ) 

sixth-order methods smaller C values are required. 

Computationally this iterative procedure has many advantages. 
The implicit inversion matrices are always strongly diagonally 
dominant. The initial spline correctors are obtained from rea- 
sonably smooth approximate solutions. In this way, it is possible 
to suppress spurious and bothersome oscillations in M ^ , etc. 

arising from the use of arbitrary initial conditions, or inaccurate 
coarse grid central-difference solutions. The final spline result 
on the same course grid may be quite reasonable. 

One of the major advantages of thio approach is improved com- 
putational efficiency. With the higher-order de ferre-V- corrector , 
it is possible to minimize the number of arithmetic operations so 
that only a marginal increase over the lower-order finite-differ- 
ence requirement results. The evaluation of the correctors are in- 
corporated in an efficient end independent subroutine. Moreover, 
the corrector curve fits need not be evaluated every iteration. 

A systematic analysis to determine the optimum corrector update 
procedure has not been rigorously considered. It has been found 
that the number of over-all iterations for convergence is not sig- 
nificantly increased if the corrector and therefore the curve fits 
are evaluated every two to three iterations. It is through this 
optimization mechanism for achieving computational efficiency that 
the present application of higher-order deferred-correction tech- 
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niqucs bocomos even more attractive. It is interesting that the 

rate of convergence of this method is comparable to that f'^'und 

in the earlier coupled spline calculations. 

The higher-order approximations of the convective terms may 

sometimes lead to an increase in aliasing error and thus non- 

(81 

linear instability. The present iterative procedure will not 
alleviate this form of instability. In such cases, additional 
considerations of smoothing and or filtering are required. 

These will be introduced, as necessary, in the present study in 
order to insure that the aliasing error remains bounded. The 
spline 4 solution for the flow in a channel with a backward step 
is one example in which this mode of instability is found. The 
introduction of smoothing in the spline curve fits, maintains the 
over-all accuracy of the numerical scheme but eliminates the non- 

/ O \ 

linear instability. This procedure will be discussed for the 
channel problem. 

4. Boundary conditions 

For the spline calculations considered here, higher-order 
boundary conditions are required for ip and w; in addition, appro- 
priate values for the corrector spline curve fits etc., 

must be specified. The various boundary conditions are as follows: 
i) At a solid surface the stream function \p is prescribed; the 
vorticity is evaluated from the no-slip condition, so that 


u> 


1,1 - ^ <*i,2-»i,l> + 


r~- {K^ 3 

6o(l+o)^ ' 


(7a) 


the portions of (7a) containing are coupled directly into the 2x2 
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algorithm, sec Part 1. The highor-ordcr spline correction in 
(7a) is treated explicitly. This is consistent with the deferred- 
correction treatment of the equations. 

ii) At a far-field boundary, the conditions Uy-*0 and <(»y-*l are 
specified. With equations (1) these conditions are of the form: 





i,N-l 






(7b) 


A similar expression obtains for As before, the K correctors 

arc treated explicitly and the ip terms implicitly. 

iii) At the inflow boundary the values of and a> are specified. 

iv) At the outflow boundary the conditions ■ 0 are im- 

posed, either directly through extrapolation or indirectly by ap- 
plying the boundary layer form of the governing equations. If 

I 

extrapolation is used, the boundary values for ip and w are coupled 
implicitly into the solution algorithm. If the boundary- layer 
equations ar'^ 'i.?. „umed at the outflow, the spline 4 deferred-correc- 
tor procedure of section 3 is used to couple the solution of these 
equations to the interior flow. This is carried out by appropriate 
modification of the algorithm coefficients. There is no upstream 
influence from the outflow boundary v.’ith the boundary layer condi- 


ti<" . 8 . 


Finally, the spline boundary conditions required for the 
curve fits of the M's, L's, etc., are obtained by satisfying the 
governing equations (3,5) at the surface; at the inflow/outflow 
the derivative conditions K, . » K^. and K„. “ , ., etc., are 

specified. These conditions on the third derivative are loss 
rigid than the direct specification of K. In many cases this leads 
to the suppression of high-frequency oscillations in the curve fits. 
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5. Examples 

The various flows examined in Part i have been reinvestigated 
with the more accurate highei-order spline 4 deferred-corrector 
strongly implicit procedure. In all the cases, solutions are ob- 
tained directly with At ■ 10^. For the flow in a channel with a 
backward facing step, the spline deferred correction solution ex- 
hibits oscillations that lead to a non-linear instability. This 
is an aliasing effect and can be eliminated by the introduction of 
filtering or smoothing of the spline curve fits in the axial di- 
rections. This effect will be detailed in a following section. 

5.1 Flow in a Driven Cavity 

This problem has been considered by many investigators using 
a variety of numerical techniques. It serves a particularly useful 
purpose here as the spline 4 solutions have already been obtained 
with an uncoupled time-dependent method and therefore the ef- 
fectiveness, stability and accuracy of the present coupled implicit 
corrector method can be assessed. Buneman's direct solver was ap- 
plied in previous calculations'^^ for the stream function and a 
predictor- corrector procedure for the vorticity equation was ap- 
plied. Significantly, the present coupled strongly implicit method 
converges considerably faster than the uncoupled direct solver/pre- 
dictor-corrector combination. An exact time comparison of the 
present procedure with the earlier spline 4 calculations is not 

possible as different computers were used. The uncoupled solutions 

(3) 

typically required about 4 minutes on a CDC 6600 computer for 
= 400 and a 17x17 grid. The coupled strongly implicit deferred- 
corrector method converges in about 75 secs on an IBM 360/65; this is 
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about 3 to 8 times slower than the 6600. Furthermore, the present 
analysis allows for a solution of the steady state form of the 
equations. The final results, for a uniform 17x17 grid, are 
identical with those presented in reference 3. The number of 
iterations for the higher-order corrector method is approximately 
the same as was required for the finite-difference solution of 
Part 1. Since the spline 4 calculations require only one-sixteenth 
the number of mesh points (in 2-dimensions) as the equally accurate 
finite-difference solutions of Part 1, the higher-order spline 
computations require much less computer time. Furthermore, the 
deferred-corrector method is relatively simple to implement and 
the spline corrector can be added as a subroutine to an existing 
finite-difference . 

5.2 Heat Transfer in a Driven Cavity 

Once converged velocity profiles are obtained, the solution of 
the energy equation 

i T^. + uT^ + vTy = ^ , (7) 

can be determined independently. The spline 4 discrete form of 
equation (7) , utilizing the iterative procedure outlined in section 
3, is the same as that for vorticity u) given by equation (4a); how- 
evej , the non-conservation form of the energy equation is used 
here . 

The thermal boundary conditions and the origin of the pres- 
cribed velocity distribution are explained in Part 1. All of the 
thermal boundary conditions are implicitly coupled into Stone's 
solution algorithm. The spline boundary conditions arc obtained by 
satisfying equation (7) at the walls of the cavity. Tlie solution 
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for a Pcclot number ■ 50 has been obtained on a 15x15 uniform 
grid. The heat transfer at the side walls is shown in figures 3,4 
where the finite difference solution of Part 1 is also reproduced 
for comparison purposes. The temperature distribution at the up- 
per and lower walls, as well as the mid-plane, is shown in fig- 
ure 2. It can be seen that the spline 4 temperature distribution 

on the upper moving wall is considerably different from that ob- 

* 

tained with the finite-difference procedure. Notably, there are 
no oscillations in the spline solutions. 

5.3 Flow in a Channel with a Backward Facing Step 

The conformal transformation that maps the step geometry into 
a straight channel, and the prescribed grid distribution are dis- 
cussed in Fig. 5. The spline 4 iterative procedure is initiated 
after 10 to 20 iterations of the central finite-difference calcu- 
lation. At the outflow boundary, the spline 4 solution to the 
i|>-(i) form of the boundary-layer equations was coupled into the 2x2 
strongly implicit algorithm for the Navier-Stokes equations. In 
order to tost the higher-order accuracy of the spline 4 method 
the fully developed flow in a straight channel was evaluated. The 
change from stcond-to fourth-order accuracy was evident. 

When the spline corrector was switched on, it was soon ob- 
served that the higher-order solution near the outflow boundary 
acquired a small oscillation. This subsequently developed into an 
instability. A number of reasons for this instability are possible. 
The most likely was the sudden transition from the fourth-order ac- 
curate representation of (uuO in the interior to the first- order 
accurate boundary layer boundary condition at the outflow. In 
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order to make this transition consistent a variety of boundary 
conditions for the spline curve fits wore examined; however, 
the instability persisted. Subsequently, systematic-numerical 
experimentation indicated that the source of this instability 
was closely related to the ip curve fit for the evaluation of 
V and the curve fit of required in the deferred corrector 

in equation (4a). Based on some previous studios' of non- 
linear instability it was apparent that the root cause was the 
' aliasing error due to the higher-order representation of the con- 
vective terms along the channel and the extremely coarse grid. 

The instability was eliminated by introducing smoothing in the 
curve fits. For example, after is evaluated, a new is ob- 
tained from the three point smoothing formula 

'^141 + + ^i ^ ‘^l^i-l 




TTTkJTT 


( 8 ) 


It can be shown that for a variable grid, formal third-order ac- 
curacy of the convective as well as diffusive terms results. 

For uniform grids the smoothing is fourth-order accurate for \p,u). 
Converged spline 4 solutions for * 1, 100 and 1000 were now ob- 
tained. The vorticity along the lower and upper walls is shown 
in figs. 6,7 where the finite-difference solution on the same grid 
is also depicted. 

With the present refined grid, the finite-difference solution 
is quite accurate. Outside of the recirculation region, tlie 
changes in the flow variables arc not large and therefore the fin- 
ite-difference approximation is also reasonable in the fully de- 
veloped flow. In the recirculation region wlierc the grid is very 
fine the spline 4 solution shows only a marginal difference. The 
peak wall vorticity is slightly larger. Also, the point of separa- 
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tion is Blightly different from that of the finite-difference 
calculations. Finally, for the fully developed region the 
spline 4 results agree with the known exact solution to three 
decimal places. For ■ 1, the flow separates below the step. 

The separation point moves towards the corner as is increased. 

SUMMARY 

A spline 4 deferred corrector strongly implicit method for 
the 4»-o) form of the Navier-Stokes equation has been investigated. 
In this iterative approach, both convective and diffusive terms 
in the governing equations are reordered such that at each level 
of the iteration a modified finite-difference form of the govern- 
ing equations is solved. The higher-order spline correction is 
therefore treated explicitly. This approach has many distinct 
advantages. For example, (i) the computer program does not in- 
volve any significant additional complexity other than would occur 
for a second-order finite-difference program. The implicit block 
size remains unchanged. ii) The additional computational effort 
is associated solely with the evaluation of the explicit spline 
curve fits. Therefore, the existing finite-difference codes can 
easily be upgraded to higher-order accuracy. iii) The over-all 
storage for the fourth-order method is considerably reduced, 
iv) Unlike the conditional stability of the method of deferred cor- 
rection (ref. 5), the present technique is unconditionally stable. 

Finally, with this coupled strongly implicit procedure, 
spline 4 "steady-state" solutions (with At = 10 have been obtained 
for a variety of flow problems. It has been found that the number 
cf iterations for the spline 4 deferred corrector procedure is com- 
parable to that required for the finite-difference solutions pre- 
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sentcd in part 1. The method of solution converges rapidly when 
compared with some of the earlier uncoupled methods applied for 
the solution of the Navier-Stokes equations. 
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